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Graphene, a hexagonal sheet of sp^-bonded carbon 
atoms, has extraordinary properties which hold immense 
promise for future nanoelectronic applications. Unfor- 
tunately, the popular preparation methods of microme- 
chanical cleavage and chemical exfoliation of graphite do 
not easily scale up for application purposes. Epitaxial 
graphene provides an attractive alternative, though there 
are many challenges, not least of which is the absence 
of an understanding of the complex atomistic assembly 
kinetics of graphene. Here, we present a simple rate the- 
ory of epitaxial graphene growth on close-packed metal 
surfaces. Based on recent low-energy electron-diffraction 
microscopy experiments (LEEM)^, our theory supposes 
that graphene islands grow predominantly by the addi- 
tion of five-atom clusters, rather than solely by the cap- 
ture of diffusing carbon atoms. With suitably chosen ki- 
netic parameters, our theory produces a time-dependent 
carbon adatom density that is in quantitative agreement 
with measured data. The temperature-dependence of 
this adatom density at the onset of nucleation leads us 
to predict that the smallest stable precursor to graphene 
growth is an immobile island composed of six five-atom 
clusters. Our findings provide a starting point for more 
detailed simulations which will yield important input to 
developing strategies for the large-scale production of epi- 
taxial graphene. 

The epitaxial growth of graphene on transition-metal 
surfaces has been studied by surface scientists for many 
yearsS.. The broader scientific community took notice 
of graphene in 2004, when two groups independently 
pointed out its unusual transport properties and poten- 
tial for microelectronic applications-i^'. Subsequent work 
has revealed additional remarkable properties, such as 
anomalies in the integer quantum Hall effect'^ and in 
quasi-particle coupling^, which are signatures of charge 
carriers that behave as massless Dirac fermions. 

The exfoliation technique of Ref. [1] is often used to 
produce samples of graphene for two-dimensional elec- 
tron gas studies. This method is quick and easy, but it 
is uncontrolled and unsuitable for scale-up to the pro- 
duction volumes that will be necessary for the most in- 
teresting and important commercial applications. An at- 
tractive alternative is to grow graphene epitaxially on 
a hexagonal substrate. Considerable progress has been 
achieved in this direction, both using silicon carbide^ 
and close-packed metals- as the substrate. Nevertheless, 
almost nothing is known about the atomistic formation 
kinetics of epitaxial graphene on any substrate. 



The present paper is motivated by particularly elegant 
LEEM experiments reported by Loginova and co-workers 
for the growth kinetics of the Gr/Ru(0001) systemic. Us- 
ing electron reflectivity to monitor the evolution of the 
density of carbon adatoms deposited either from a heated 
carbon rod or from an ethylene source, it was discovered 
that the step edges of two-dimensional (2D) graphene is- 
lands advance with a velocity v that varies as the fifth 
power of the ratio of the carbon adatom density n to its 
equilibrium value noq: 
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This led the authors to suggest that the islands grow 
by the attachment of five-atom clusters, rather than by 
the usual mechanism of single adatom attachment. More 
recent work by the same group extends this conclusion 
to the Gr/Ir(lll) system-'^. 

In this paper, we develop a rate theory for the epi- 
taxial growth of graphene on metal surfaces. Rate equa- 
tions have been used for some time for modelling epi- 
taxial systemsii and for establishing general principles 
such as the existence of scaling regimes for concentra- 
tions of surface species and scaling forms for island-size 
distributions^^. Scanning tunneling microscopy images, 
both as "snapshots" of quenched surfaces and in situ 
scansi^, are usually interpreted in terms of kinetic Monte 
Carlo simulations because of the wealth of spatial infor- 
mation contained in such images^i^. Nevertheless, rate 
equations continue to provide the conceptual framework 
for interpreting essentially all growth scenarios. The 
amount and quality of time-resolved data in the present 
case is unprecedented and, as we will show, the temper- 
ature dependence of the carbon adatom density at the 
onset of nucleation and at equilibrium encode key infor- 
mation about the atomistic processes of graphene forma- 
tion on Ru(0001)i that would be unavailable without a 
theoretical analysis. 

A simple model based on the experiments of Loginova 
et al}^ assumes that carbon atoms arrive at the surface 
with a flux F and migrate across the surface with diffu- 
sion constant D. Adatoms can attach to the perimeter of 
2D islands, but this is not the predominant mechanism 
of island growth. Instead, clusters composed of i = 5 
carbon atoms form when i adatoms collide. These clus- 
ters migrate with surface diffusion constant D' and, when 
j — (!) such clusters collide, an island of size i x j is formed. 
Islands grow mainly by the attachment of i-atom clusters. 
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The value i = 5 for the size of mobile clusters is derived 
from experiment, as noted above, while the value j = 6 
is obtained from a fit to the temperature dependence of 
the carbon adatom density at the onset of nucleation. 
The latter is a new and striking result, and suggests that 
the collision of fewer than six clusters leads to the forma- 
tion of transient, possibly mobile, species whose struc- 
ture bears little resemblance to graphene. Even islands 
formed from six mobile clusters are unlikely to immedi- 
ately transform into graphene, but instead pass through 
a series of intermediate structures. This general scenario 
is consistent with tight-binding, grand canonical Monte 
Carlo simulations which produce small carbon chains as 
a precursor to graphene formation on Ni(lll}^. 

The rate equations for these kinetic steps are expressed 
in terms of the homogeneous densities of carbon adatoms 
n, five-atom clusters c, and islands N as 



= F- iDn' + iKc - DnN + K'N . 



dn 
'dt 
dc 

— = Dn' -Kc- D'cN - iD'c' 
dt 
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= D'c^. 
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In these equations, K is the cluster dissolution rate and 
K' is the detachment rate of adatoms from islands. For 
simplicity, we have omitted refinements to these equa- 
tions, such as capture numbers, that make them more 
realistioi^^, but also more numerically cumbersomel^. 
Thus, our proposed rate equations are sufficiently mini- 
mal to enable a straightforward optimization of param- 
eters, while retaining enough physical content to permit 
a detailed interrogation of experimental data. We will 
demonstrate below that this set of equations is capable 
of providing a quantitative account of all of the measured 
data reported in Ref. T]. Moreover, our analysis will re- 
veal the presence of new kinetic mechanisms for graphene 
formation on metal surfaces. 

Each of the rates parameters D, D' , K, and K' is taken 
to have the Arrhenius form i'oe~^^ with the prefactor 
h'o = 2kBT/h assumed to be the same for all processes. 
Otherwise, ks is Boltzmann's constant, T is the absolute 
temperature of the substrate, h is Planck's constant, /3 = 
l/{kBT), and E is the energy barrier to the process. For 
the temperature range of interest (790 K-1050 K), i/q ^ 
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We have set the lattice constant a = 1, but 



the required factors of a can be reinstated into (HJ-Q 
by straightforward dimensional analysis. 

There are four energy barriers that must be deter- 
mined, as well as the value of j. However, this is not 
an unconstrained optimization, as there are restrictions 
which limit the ranges of these parameters. For exam- 
ple, the equilibrium carbon adatom concentration ricq 
[Fig. (Tfa)], determined from the time-independent solu- 
tion to © with -F = 0, yields 



(5) 



As rieq is found to be an increasing function of 
temperature^, we must choose Ek' > Ejj. Indeed, by 
fitting this form of ricq to the corresponding experimen- 
tal values, we find that Ek' — Ed — 0.35 eV. In obtaining 
(O, we have used the fact that n ^ c ^ N , which results 
from the unusually large values of i and j. This can be 
checked afterward for self-consistency (Fig. [ij. 

Further insight into the growth parameters is obtained 
by examining the stationary concentration in the pres- 
ence of a non-zero flux. Again invoking the ordering 
n ^ c ^ TV, we find. 
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where — Ft is the coverage, 
perature-dependence of rig is 



At fixed 9, the tem- 
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Ed' -jEk + {j + 1)Ed 



(7) 



As the behavior with temperature of Us is qualitatively 
similar to that of Tinuc, which is observed to be an in- 
creasing function of temperature^, we must require that 



Ed' - jEk + {] + 1)Ed < . 



(8) 



Note that this condition involves j but not i. For our 
energy parameters, we must choose j > 2. Beginning 
with five parameters (-B_d, -Eds Ek, Ek', and j), the 
restrictions in ([5]) and ([8]) reduce our optimization to an 
effective three-parameter fit. 

The solutions of the rate equations for the evolutions 
of the densities of the carbon adatoms, the five-atom 
clusters, and the immobile islands are shown in Fig. [T] 
with the optimized parameters in Table HI These solu- 
tions are most sensitive to the energy barriers for D, K , 
and K', with variations of ±0.01 eV producing notice- 
able changes in the solutions. In contrast, variations by 
as much as ±0.05 eV of the diffusion barrier D' for the 
five-atom clusters produce minimal changes, especially 
for the adatom solution. This is due, in part, to the 
fact that the processes that most influence the adatom 
density are the formation and dissolution of the clusters. 
The density of islands and their sizes have a smaller ef- 
fect because the attachment and detachment processes 
between adatoms and islands occur much less frequently, 
as the island step velocity in ([T]) indicates. 

Most apparent from the solutions in Fig. [T] is how 
qualitatively different the adatom and island profiles are 
from those found in the "usual" epitaxial growth scenario 
where only adatoms are mobile and there is a critical size 
for island formation^^-^^. The presence of mobile clusters 
composed of five atoms and the large size (j — 6) of the 
smallest immobile islands formed from five-atom clusters 
are responsible for this behavior. The clusters and is- 
lands both show sudden increases in their densities due 
to bursts of aggregation when there is a sufficient adatom 
density to first form five-atom clusters. The onset of is- 
land formation is further delayed until there are enough 
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five-atom clusters so that a six-cluster collision becomes 
likely. The islands are formed over a short time period 
after which the island density sasturates with essentially 
no additional nucleation. 

Figure [2] shows the solution for the carbon adatom den- 
sity in Fig. [ija) overlayed with the experimental data in 
Fig. 12 of Ref. The theoretical solution reproduces 
all of the main features of the experimental data, though 
there are small differences in some of the details. The rate 
equations account for the position of the nucleation peak 
at approximately 150 s after the flux is introduced into 
the system, but underestimate the experimental value of 
JT-nuci by ~ 10%. The subsequent stationary regime is 
also accounted for by our theory, but the solution overes- 
timates the steady experimental value by a few per cent. 
After the cessation of the flux, the theoretical and exper- 
imental equilibration curves are indistinguishable. 

The most likely cause of the discrepancies between the 
theoretical and experimental carbon adatom evolutions is 
the omission of any spatial information in the rate equa- 
tions. By using spatially averaged densities for all surface 
species, the local spatial arrangements of five-atom clus- 
ters and islands, which affects adatom attachment and 
detachment through the notion of "capture zones" 
and any preferential sites for cluster or island formation, 
are excluded. These would first be apparent at the on- 
set of nucleation, where the spatial arrangement of is- 
lands is critical to establishing their growth ratesi^. The 
subsequent discrepancy in steady-state regime is then a 
consequence of this underestimate of the nucleation rate. 

On the other hand, the rate equations provide an ex- 
cellent account of the equilibration of the adatom con- 
centration after the cessation of the incident mass flux. 
The comparative immunity of the equilibration regime 
to spatial correlations between islands suggests that n^q 
is dominated by the detachment of adatoms from an is- 
land followed by re-attachment to the same island. An 
immediate consequence of this is that the distribution 
of island sizes would not change over the time scale of 
equilibration. This can be tested experimentally. 



Figure[3]compares the temperature-dependence of Hnuc 
and TT-eq [Fig.[IJa)] obtained from the rate equations with 
those observed in experiments^. As Fig. [5] indicates, 
is expected to be more accurate than nnuc, although 
both quantities provide good accounts of the experimen- 
tal data. The computed equilibrium concentrations from 
the full solution are indistinguishable from that in ([S]), 
which shows that, while single adatom attachment and 
detachment do not have a significant direct impact on 
island kinetics during growth, they do have experimen- 
tally observable consequences during equilibration. On 
the other hand, the temperature-dependence of rij^uc em- 
bodies more information about kinetics during growth, 
as Eq. d?]) indicates, especially with regard to the num- 
ber of five-atom clusters required to form an incipient 
graphene island. The choice j = 6 yields the best fit 
to the data in Fig. |3l This comparison, more than any 
other reported here, indicates the value of our rate equa- 
tion analysis. An increasing rinuc with temperature is not 
at all typical of epitaxial systems, and its consequences 
would have been difficult, if not impossible, to ascertain 
without a theoretical model that abstracts all but the 
most essential kinetic steps. 

Methods Rate equations are typically formulated phe- 
nomenologically for a particular system with basic guid- 
ance from experiments and fundamental calculations. In 
this case, we began with processes involving the diffusing 
carbon adatoms and five-atom clusters. Additional pro- 
cesses and all parameters were identified through com- 
parison with the available experimental data. The math- 
ematical structure of rate equations for epitaxial systems 
is a system of coupled autonomous nonlinear ordinary dif- 
ferential equations that often have a property known as 
"stiffness," which arises from the disparity of rates in the 
equations. This means that their solutions become unsta- 
ble in certain parameter regimes. Our calculations were 
carried out with Mathematical^, which has a solver 
that can accommodate such cases. 
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FIG. 1: Solution of rate equations. Coverage dependence of (a) the carbon adatom density n, (b) the five-atom mobile 
cluster density c, and (c) the total immobile island density obtained from the rate equations ©-(SI) at T = 1020 K. The 
adatom densities at the onset of nucleation nnuc and at equilibrium n^q are indicated in (a). The kinetic parameters are 
Ed = 0.92 eV, Eo' = 0.87 eV, Ek = 1.72 eV, Ek' = 1.27 eV, and j = 6. 
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FIG. 2: Rate theory and LEEM measurements for time-dependent carbon adatom coverage. The evolution of the 
carbon adatom density in Fig. [2a) (black curve) compared with the corresponding LEEM data in Ref [l|] (red points) 
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FIG. 3: Rate theory and LEEM measurements for adatom coverage at onset of nucleation and at equilibrium. 

The temperature dependence of the concentration of carbon adatoms nnuc required to nucleate graphene (red symbols) and 
the carbon adatom concentration ricq at equilibrium (blue symbols) compared with the corresponding quantities obtained from 
rate equations (black symbols) as indicated in Fig. [T] Experiments^ were carried out with carbon vapor (filled red and blue 
symbols) and ethylene (open red and blue symbols). 
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TABLE I: The optimized energy barriers (in eV) in the Arrhenius forms of D, D' , K, and K' in Q-®- The accuracy of each 
barrier refers to the sensitivity of the solutions to variations of that barrier. 



Rate Barrier Accuracy 

D 0.92 0.01 

D' 0.87 0.05 

K 1.72 0.01 

K' 1.27 0.01 



